!-------------------------------------------------------------------------------------
! subroutine UpWind2ndOrder
!-------------------------------------------------------------------------------------

subroutine UpWind2ndOrder(  u,  dt, dx, a, jx, uLeft,  uRight )
implicit none
integer jx, i
real u(jx),  dt, dx, a, uLeft,  uRight
real sigma 
! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
! du/dt = -a du/dx
! For a >0

! sigma = a * dt / dx ,  sigma < 2 for stability

! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

!            all u's on RHS are at time step "n"

sigma = a * dt / dx

!write(*,*) 'sigma=', sigma,'   jx=',jx  ; pause 'pause 01 in subroutine'
!write(*,*)'u=', u     ; pause 'pause 02 in subroutine'

if( sigma > 0.0 ) then  
       do i=jx,3,-1
         u(i) = u(i) - (sigma/2)*( 3*u(i) - 4*u(i-1) + u(i-2) ) + (sigma*sigma / 2 )*( u(i) - 2*u(i-1) + u(i-2) )
       enddo
       u(2) = u(2) - sigma*( u(2) - u(1) )!  first order for u(2)
       u(1) = uLeft
else  
       do i=1, jx-2
         u(i) = u(i) + (sigma/2)*( 3*u(i) - 4*u(i+1) + u(i+2) ) + (sigma*sigma / 2 )*( u(i) - 2*u(i+1) + u(i+2) )
       enddo
       u(jx-1) = u(jx-1) - sigma*( u(jx) - u(jx-1) )!  first order for u(2)
       u(jx) = uRight
endif
         


return;end


!-------------------------------------------------------------------------------------
! subroutine UpWind2ndOrder_Variable_a
!-------------------------------------------------------------------------------------

subroutine UpWind2ndOrder_Variable_a( u, dt, dx, aVar, jx, uLeft,  uRight, iPeriodic )
! as per communication from B. Wendroff on Jan 10, 2011, for:  du/dt =  - (d/dx)[ aVar(x) u ]
! see Word doc:   Advection in non uniform flow
implicit none
integer jx, i, iPeriodic
real u(jx), dt, dx, aVar(jx), uLeft,  uRight
real sigma
real    u0,    uM1,    uJXp1,    uJXp2
real aVar0, aVarM1, aVarJXp1, aVarJXp2
! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
! du/dt = -aVar du/dx for aVAr constant
! sigma = aVar_max * dt / dx ,  sigma < 2 for stability if aVar is uniform

! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

!            all u's on RHS are at time step "n"

  sigma =     dt / dx

     u0 =    u(jx);    uM1 =    u(jx-1);    uJXp1 =    u(1);    uJXp2 =    u(2)  ! needed for periodic boundary conditions
  aVar0 = aVar(jx); aVarM1 = aVar(jx-1); aVarJXp1 = aVar(1); aVarJXp2 = aVar(2)  ! needed for periodic boundary conditions

if( aVar(1) > 0.0 ) then   !  this implicitly assumes that aVar is mostly positive, or else, mostly negative.
       do i=jx,3,-1
         u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i-1)*u(i-1) + aVar(i-2)*u(i-2) ) + (sigma*sigma / 2 ) * &
( 0.5*( aVar(i-1)+aVar(i) )*( aVar(i)*u(i)-aVar(i-1)*u(i-1) ) - 0.5*( aVar(i-2)+aVar(i-1) )*( aVar(i-1)*u(i-1)-aVar(i-2)*u(i-2) )  )
               
!  {          aVar(i-1/2)  * [  aVar(i)u(i)-aVar(i-1)u(i-1) ]  -         aVar(i-3/2)      *  [ aVar(i-1)u(i-1)-aVar(i-2)u(i-2) ]   }
       enddo
       if ( iPeriodic .NE. 1) then
          u(2) = u(2) - sigma*( aVar(2)*u(2) - aVar(1)*u(1) )!  first order for u(2)
          u(1) = uLeft
       else ! periodic BC
          i = 2
          u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i-1)*u(i-1) + aVar0*u0 ) + (sigma*sigma / 2 ) * &
    ( 0.5*( aVar(i-1)+aVar(i) )*( aVar(i)*u(i)-aVar(i-1)*u(i-1) ) - 0.5*( aVar0+aVar(i-1) )*( aVar(i-1)*u(i-1)-aVar0*u0 )  )
          i=1
          u(i) = u(i) - (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar0*u0 + aVarM1*uM1 ) + (sigma*sigma / 2 ) * &
          ( 0.5*( aVar0+aVar(i) )*( aVar(i)*u(i)-aVar0*u0 ) - 0.5*( aVarM1+aVar0 )*( aVar0*u0-aVarM1*uM1 )  )          
       endif
else !  aVar(1) < 0 
       do i=1, jx-2
         u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i+1)*u(i+1) + aVar(i+2)*u(i+2) ) + (sigma*sigma / 2 ) * &
( 0.5*( aVar(i+1)+aVar(i) )*( aVar(i)*u(i)-aVar(i+1)*u(i+1) ) - 0.5*( aVar(i+2)+aVar(i+1) )*( aVar(i+1)*u(i+1)-aVar(i+2)*u(i+2) )  )
       enddo
       if( iPeriodic .NE. 1) then
         u(jx-1) = u(jx-1) - sigma*( aVar(jx)*u(jx) - aVar(jx-1)*u(jx-1) )!  first order for u(jx-1)
         u(jx) = uRight
       else ! periodic BC
         i = jx - 1
         u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVar(i+1)*u(i+1) + aVarJXp1*uJXp1 ) + (sigma*sigma / 2 ) * &
   ( 0.5*( aVar(i+1)+aVar(i) )*( aVar(i)*u(i)-aVar(i+1)*u(i+1) ) - 0.5*( aVarJXp1+aVar(i+1) )*( aVar(i+1)*u(i+1)-aVarJXp1*uJXp1 )  )
         i = jx
         u(i) = u(i) + (sigma/2) * ( 3*aVar(i)*u(i) - 4*aVarJXp1*uJXp1 + aVarJXp2*uJXp2 ) + (sigma*sigma / 2 ) * &
         ( 0.5*( aVarJXp1+aVar(i) )*( aVar(i)*u(i)-aVarJXp1*uJXp1 ) - 0.5*( aVarJXp2+aVarJXp1 )*( aVarJXp1*uJXp1-aVarJXp2*uJXp2 )  )
       endif
endif
         
return;end


!-------------------------------------------------------------------------------------
! subroutine SubUpwind_NonCon_Form
!-------------------------------------------------------------------------------------


subroutine SubUpwind_NonCon_Form( u, dt, dz, aVar, jz, uBottom,  uTop, iPeriodic, iUpWind, uTemp )
! as per communication from B. Wendroff on Jan 10, 2011, for:  du/dt =  - (d/dz)[ aVar(x) u ]
! see Word doc:   Advection in non uniform flow
! this version is similar to that in "linear lib" folder but source has been dropped
implicit none
integer jz, i, iPeriodic, iUpWind
real u(jz), dt, dz, aVar(jz), uBottom,  uTop, uTemp(jz)
real sigma
real    u0,    uM1,    uJXp1,    uJXp2
real aVar0, aVarM1, aVarJXp1, aVarJXp2
! Vol. 2 page 502 of Numerical Computation of Internal and External Flows by Hirsch
! du/dt = -aVar du/dz for aVAr constant
! sigma = aVar_max * dt / dz ,  sigma < 2 for stability if aVar is uniform

! u(i,n+1) = u(i) - (sigma/2)*[ 3*u(i) - 4*u(i-1) + u(i-2) ] + (sigma*sigma / 2 )*[ u(i) - 2*u(i-1) + u(i-2) ]

!            all u's on RHS are at time step "n"

  sigma =     dt / dz

     u0 =    u(jz);    uM1 =    u(jz-1);    uJXp1 =    u(1);    uJXp2 =    u(2)  ! needed for periodic boundary conditions
  aVar0 = aVar(jz); aVarM1 = aVar(jz-1); aVarJXp1 = aVar(1); aVarJXp2 = aVar(2)  ! needed for periodic boundary conditions
  
IF( iUpWind==1 ) THEN

if( aVar(1) > 0.0 ) then   !  this implicitly assumes that aVar is mostly positive, or else, mostly negative.
       do i=jz,3,-1
         u(i) = u(i)      -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i-1) + u(i-2) )  &
              + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i-1) + u(i-2) )  &
              +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i-1) + u(i-2) )  * ( 3*aVar(i) - 4*aVar(i-1) + aVar(i-2) )             
       enddo
       if ( iPeriodic .NE. 1) then
          i = 2
          u(2) = u(2) - sigma * aVar(2) * ( u(2) - u(1) )!  first order for u(2)
          u(1) = uBottom
       else ! periodic BC
          i = 2
         u(i) = u(i)    -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i-1) + u0 )  &
            + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i-1) + u0 )  &
            +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i-1) + u0 )  * ( 3*aVar(i) - 4*aVar(i-1) + aVar0 )  
          i=1
         u(i) = u(i)    -   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u0 + uM1 )  &
            + (sigma*sigma / 2 ) * aVar(i)**2  * (  u(i)  - 2*u0 + uM1 )  &
            +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u0 + uM1 )  * ( 3*aVar(i) - 4*aVar0 + aVarM1 )   
       endif
else !  aVar(1) < 0 
       do i=1, jz-2
         u(i) = u(i) +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i+1) + u(i+2) )  &
         + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) - 2*u(i+1) + u(i+2) )  &
         +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i+1) + u(i+2) )  * ( 3*aVar(i) - 4*aVar(i+1) + aVar(i+2) ) 
       enddo
       if( iPeriodic .NE. 1) then
         i = jz - 1
         u(jz-1) = u(jz-1)-sigma * aVar(jz-1)*( u(jz) - u(jz-1) )!  first order for u(jz-1)
         u(jz) = uTop
       else ! periodic BC
         i = jz - 1
         u(i) = u(i)    +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*u(i+1) + uJXp1 )  &
            + (sigma*sigma / 2 ) * aVar(i)**2  * (   u(i) -2*u(i+1)  + uJXp1 )  &
            +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*u(i+1) + uJXp1 )  * ( 3*aVar(i) - 4*aVar(i+1) + aVarJXp1 )  
         i = jz
         u(i) = u(i) +   (sigma/2) *aVar(i) * ( 3*u(i) - 4*uJXp1 + uJXp2 )  &
         + (sigma*sigma / 2 ) * aVar(i)**2  * ( u(i)   - 2*uJXp1 + uJXp2 )  &
         +     (sigma*sigma / 8 ) * aVar(i) * ( 3*u(i) - 4*uJXp1 + uJXp2 )  * ( 3*aVar(i) - 4*aVarJXp1 + aVarJXp2 ) 
       endif
endif

ELSE  ! else use BW recommended central difference  ( modified by Source term )


do i = 2, jz-1
       uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( u(i+1) - u(i-1) ) &
                        + (sigma*sigma / 4 ) * aVar(i) * (  (aVar(i)+aVar(i+1))*(u(i+1)-u(i)) - (aVar(i)+aVar(i-1))*(u(i)-u(i-1))  )
enddo
       if ( iPeriodic .NE. 1) then
          if( aVar(1) > 0 )then
             uTemp(1) = uBottom
             i=jz
             uTemp(i) = u(i)  - sigma * aVar(i) * ( u(i) - u(i-1) )
          else
             uTemp(jz) = uTop
             i = 1
             utemp(i) = u(i)  - sigma * aVar(i) * ( u(i+1) - u(i) )
          endif
       else ! periodic
       i=1
       uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( u(i+1) - u0 ) &
                        + (sigma*sigma / 4 ) * aVar(i) * (  (aVar(i)+aVar(i+1))*(u(i+1)-u(i)) - (aVar(i)+aVar0)*(u(i)-u0)  )
       i=jz
       uTemp(i) = u(i)  - (sigma/2) *aVar(i) * ( uJXp1 - u(i-1) ) &
                        + (sigma*sigma / 4 ) * aVar(i) * (  (aVar(i)+aVarJXp1)*(uJXp1-u(i)) - (aVar(i)+aVar(i-1))*(u(i)-u(i-1))  )
       endif   !  end of periodicity branch

       u(:) = uTemp(:)

ENDIF  !  end of upwind versus centered choice
         
return;end


!-------------------------------------------------------------------------------------
! subroutine NonLinearAdvection
!-------------------------------------------------------------------------------------

subroutine NonLinearAdvection( u, dt, dx, jx, uLeft,  uRight, iPeriodic )
! see Word doc:   Advection in non uniform flow.doc
!  /Users/harveyrose/Documents/MyLib/Advection in non uniform flow.doc
!  advances one time step, dt,  du/dt = - u du/dx
implicit none
integer jx, i, iPeriodic
real u(jx), dt, dx,  uLeft,  uRight
real sigma
real    u0,    uM1,    uJXp1,    uJXp2

       u0 =    u(jx);    uM1 =    u(jx-1);    uJXp1 =    u(1);    uJXp2 =    u(2)  ! needed for periodic boundary conditions

  sigma =     dt / dx
  
  IF( u(1) > 0 ) THEN
       do i=jx,3,-1
         u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i-1)**2 + u(i-2)**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( u(i-2)**3 + u(i)**3 - 2 * u(i-1)**3 )
       enddo
       if ( iPeriodic .NE. 1) then
         u(2) = u(2) - 0.5D0 * sigma * (  u(2)**2 - u(1)**2 )
         u(1) = uLeft
       else                        ! periodic and u > 0
         i=2
         u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i-1)**2 + u0**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( u0**3 + u(i)**3 - 2 * u(i-1)**3 )
         i=1
         u(i) = u(i) - 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u0**2 + uM1**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( uM1**3 + u(i)**3 - 2 * u0**3 )
         
       endif
  ELSE                             ! u < 0 
       do i=1, jx-2
         u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i+1)**2 + u(i+2)**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( u(i+2)**3 + u(i)**3 - 2 * u(i+1)**3 )
       enddo
       if ( iPeriodic .NE. 1) then
         u(jx-1) = u(jx-1) - 0.5D0 * sigma * (  u(jx)**2 - u(jx-1)**2 )
         u(jx) = uRight
       else                        ! periodic and u < 0
         i = jx-1
         u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * u(i+1)**2 + uJXp1**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( uJXp1**3 + u(i)**3 - 2 * u(i+1)**3 )
         i = jx
         u(i) = u(i) + 0.25E0 * sigma * ( 3 * u(i)**2 - 4 * uJXp1**2 + uJXp2**2 ) + &
                (1.E0 / 6.E0) * sigma * sigma * ( uJXp2**3 + u(i)**3 - 2 * uJXp1**3 )
       endif
  ENDIF
  
  return
  end

!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------

